Cox regression modeling of survival after chemotherapy for colon cancer
I. Introduction
The Cox Proportional Hazards (PH) model is one of the most widely used regression models to study survival analysis, also known as time-to-event analysis. The model investigates the association between survival time and one or more predictors. The Cox PH model is based on the hazard function, which is defined as the probability that an individual will experience the event at a given time, provided that the event has not occurred yet. The outcome variable is the survival time, which is measured from a defined time until the occurrence of an event, such as death, or the end of the study period [1]. For each subject, there is time to event (T), i.e., the time from a defined time until the event occurs, or censoring time (C), i.e., the time the subject drops out of the study or the study ends.
The Cox PH model relies on two assumptions: 1) random censoring, and 2) the proportional hazard assumption [2]. The random censoring assumption is met if patients who are censored are a random sample from the study population. There is no statistical test for this assumption, and it can be achieved through rigorous experimental design. The proportional hazard assumption states that the hazard function (hazard ratio) for two groups (e.g., experimental arm and standard arm) should remain proportional, i.e., the hazard ratio is constant over time. However, this assumption is violated when the effect of a covariate on the outcome is not constant over the follow-up period, which is common in biochemical research. When the proportional hazard assumption is not met, several methods can be applied to address this issue, including stratification of the model by variables that violate the assumption. This allows the baseline hazard to vary across strata [3]. Another approach is to include interactions between covariates and time, allowing the hazard ratios to change over time [4]. While the baseline hazards are allowed to differ, the effects of other covariates are assumed to be the same across strata. One limitation with this approach is that the effect of the stratification variable itself cannot be tested because it is not included as a covariate in the model. Other statistical methods to account for non-proportional hazards have been described and are discussed elsewhere [5].
One of the properties of the Cox PH model contributing to its popularity is that the baseline hazard, i.e., \(h_0\), is an unspecified field, which makes it a semi-parametric model. It is robust and can provide reliable estimates without specifying the baseline hazard function. Even though the baseline hazard is not known, it is possible to estimate the coefficients in the exponential part of the equation. Hence, the effect of variables included in the model can be estimated [6]. Another feature of the Cox PH model is its interpretability. The effect of covariates on the hazard is represented by hazard ratios, i.e., the relative likelihood of an event happening in the experimental group with respect to the standard group.
The Cox PH model is widely used across various fields, including medical research, epidemiology, business, engineering, and social sciences [6]. It is commonly used to determine survival rates among cancer patients with different subtypes of cancer [7], or between those treated with different treatments [8]. The Cox PH model has also been used effectively in non-medical research, for example, in the prediction of time to loan defaults [9], customer time to churn [10], product lifespan and failure time analysis in engineering [11], and for the identification of factors that influence the time it takes to achieve certain rates of success, such as vaccination rates, in public health [12].
The current analysis evaluated data from one of the first successful trials of adjuvant chemotherapy for colon cancer. The colon cancer patients in the study had their primary treatment of surgery and the objective was to test whether treatment with either Levamisole or Levamisole in combination with 5-FU chemotherapeutic agents improves survival. Levamisole is a low-toxicity compound previously used to treat worm infestations in animals; 5-FU is a moderately toxic chemotherapy agent. There are two records per person, one for recurrence and one for death. The purpose of this project is to compare survival between the untreated (Obs) group vs those treated with amisole (Lev), or amisole + 5-FU.
II. Methods
The Cox proportional hazards model was used to model the relationship between survival time and different colon cancer treatments. In particular, the survival time between the untreated group (observation) and those treated with amisole (Lev) or amisole + 5-FU was compared. The Cox regression model was chosen for this study because it is well-suited for studying the association between survival time of patients and predictors, and it allows estimating the risk or hazard of death increased or decreased due to each treatment relative to no treatment. The time (in days) until the event, i.e., death, will be modeled as a function of treatment and other variables, including age, sex, and various tumor characteristics. Significant predictors were included in the final model.
IIa. Data Source
The colon cancer cancer data set is a built-in data set in the Survival R package [13] [14]. The dataset is closest to the original study published in 1991 [15]. The Data set includes 929 subjects with stage B or C colon cancer who were randomized to three treatment groups, Observation, Levamisole (Lev), or Levamisole + Fluorouracil (5-FU). The data set includes survival and recurrence information on 929 individuals, but for the current analysis, it was filtered to include only the rows where death was the outcome variable. The time to death is given in days. The median follow-up period was 1,976 days (5.4 years), (range 23-3,329 days). The dataset includes various patient characteristics, such as demographics and tumor details, as well as the duration from surgery to trial registration, categorized as either short or long.
Column names:
| id: | id |
| study: | 1 for all patients |
| rx: | Treatment - Obs(ervation), Lev(amisole), Lev(amisole)+5-FU |
| sex: | 1=male |
| age: | in years |
| obstruct: | obstruction of colon by tumour |
| perfor: | perforation of colon |
| adhere: | adherence to nearby organs |
| nodes: | number of lymph nodes with detectable cancer |
| time: | days until event or censoring |
| status: | censoring status |
| differ: | differentiation of tumour (1=well, 2=moderate, 3=poor) |
| extent: | Extent of local spread (1=submucosa, 2=muscle, 3=serosa, 4=contiguous structures) |
| surg: | time from surgery to registration (0=short, 1=long) |
| node4: | more than 4 positive lymph nodes |
| etype: | event type: 1=recurrence,2=death |
IIb. Exploratory data analysis
The data was evaluated for missing values, duplicate entry, and outliers. The Survminer [16] package was used to plot the Kaplan-Meier curve to visualize the survival probability over time for each of the categorical variables. Multicolinearity was tested using Variant Inflation Factor (VIF) calculated using MASS package [17]. Multicolinearity was tested using Variant Inflation Factor (VIF) calculated using MASS package [17].
IIc. Statistical analysis
The R statistical software version 4.3.2 [18] was used for all analysis. The Survival package was used to construct the Cox regression model [13] [14].
Cox regression model is based on the hazard function \(h_x(t)\) with covariates at time t given by:
\(h_x(t)=h_0(t)\exp(\beta_1x_1 +\beta_2x_2 + \dots + \beta_p x_p)\)
Where:
\(h_x(t)\) is the hazard function
\(h_0(t)\) is the baseline hazard function
\(\beta_1x_1 + \beta_2x_2 + \dots +\beta_p x_p\) represent the linear combination of covariates and their coefficient
The hazards for the observation vs. amisole (Lev), or amisole + 5-FU group with covariate values x1 and x2 are given by: \(hx_1(t)=h_0(t)\exp(\beta_1x_1)\) and \(hx_2(t)=h_0(t)\exp(\beta_2x_2)\), respectively
The hazard ratio is expressed as: HR = \(hx_2(t)\) / \(hx_1(t)\) = \(\exp[\beta(x_2-x_1)]\)
The R MASS package was used for Step-wise variable selection, using “both” forward and backward variable selection [17]. For Step-wise selection, stepAIC() function uses AIC (Akaike Information Criterion) as the measure to add or remove predictors from the model.
IId. Model evaluation
Diagnostic for proportional hazard assumption
The Schoenfeld residual plot was constructed to test Cox proportional hazards assumption. When the proportional hazards assumpiton was not met for any of the covariates, stratification approach was used.
K-fold cross-validation
Model performance was evaluated using 5-fold cross-validation using the boot package [19] [20].
Model selection
The Concordance (c-index), AIC, and BIC of corresponding models were compared to select the best fit model. The model with the lowest AIC and BIC values and highest concordance was selected as the final model.
III. Analysis and Results
IIIa. Data cleaning and feature engineering
The data set contained missing values in number of nodes and differentiation columns. Missing values in differentiation and node columns were replaced with with mode and median, respectively. Numbers representing different categories in differentiation, local spread and time to surgery columns were replaced with the descriptive names. Evaluation of nodes and node4 showed that samples with greater than 4 lymph nodes have less than 5 count in nodes column, so the two columns are not consistent. Therefore, nodes column will not be used for further analysis.
IIIb. Exploratory data analysis
Figure 1: Exploratory data analysis
Exploratory data visualization shows (Figure 1.) The study participants were equally randomized in the three treatment groups, with about 33% of individuals in observation, Levamisole (Lev) or Lev + 5-FU groups. Fifty-two percent of the participants were male. The dataset is balanced with respect to the outcome (status), with 51% censored and 49% died bevore the end of follow-up period.
With respect to tumor characteristics, obstruction of colon by tumor (obstruct), perforation of colon (perfor), tumor adherence to nearby organs (adhere) was observed in 19%, 3% and 15% of the patients, respectively. About 27% of the participants have more than 4 tumor positive lymph nodes (node 4) and long time from surgery to registration (surg). Overall, most of the categories of tumor characteristics are well represented, with about 15% of the participants with those characteristics, except for perforation, differentiation and local spread.
Age is normally distributed. Evaluation of time to death or censoring showed a bimodal distribution. Stratifying the time variable by outcome (status), revealed that the majority of the patients died within 3 years of the study period.
Table 1: Patient characteristics
| Observation (%) | Amisole (%) | Amisole + 5-FU (%) | ||
|---|---|---|---|---|
| N=315 | N=310 | N=304 | ||
| Demographics | ||||
| Male | 166 (52.3) | 177 (57.1) | 141 | |
| Median age (years) [IQR] | 60 [53,68] | 61 [53,69] | 61 [52,70] | |
| Cancer characteristics | ||||
| Colon obstruction | 63 (20.0) | 63 (20.3) | 54 (17.8) | |
| Colon perforation | 9 (2.9) | 10 (3.2) | 8 (2.6) | |
| Adherence to nearby organs | 47 (14.9) | 49 (15.8) | 39 (12.8) | |
| Differentiation of tumor | ||||
| Well | 27 (8.6) | 37 (11.9) | 29 (9.5) | |
| Moderate | 236 (74.9) | 229 (73.9) | 221 (72.7) | |
| Poor | 52 (16.5) | 44 (14.2) | 54 (17.8) | |
| Extent of local spread | ||||
| Contiguous | 20 (6.3) | 12 (3.9) | 11 (3.6) | |
| Muscle | 38 (12.1) | 36 (11.6) | 32 (10.5) | |
| Serosa | 249 (79.0) | 259 (83.5) | 251 (82.6) | |
| Submucosa | 8 (2.5) | 3 (1.0) | 10 (3.3) | |
| More than 4 lymph nodes with cancer | Yes | 87 (27.6) | 89 (28.7) | 79 (26.0) |
| Short time from surgery to registration (%) | Yes | 91 (28.9) | 80 (25.8) | 76 (25.0) |
As shown in Table 1, the baseline patient characteristics were very similar between the three treatment groups, indicating suggesting that and differences in survival time between the three groups is not likely attributed to differences at the start of the study.
IIIc. Survival curves
Figure 2: Censoring (left) or survival (right) time of patients
The median survival time of patients who were censored (left) or died (right) was 2331 and 775 days, respectively. This results shows that patients who were not censored died early during the follow-up period.
Stratified survival curves
#| echo: false
#| message: false
#| warning: false
#| include: true
f1 <- survfit(Surv(time, status) ~ adhere, data = df)
f2 <- survfit(Surv(time, status) ~ rx, data = df)
f3 <- survfit(Surv(time, status) ~ sex, data = df)
f4 <- survfit(Surv(time, status) ~ obstruct, data = df)
f5 <- survfit(Surv(time, status) ~ perfor, data = df)
f6 <- survfit(Surv(time, status) ~ node4, data = df)
f7 <- survfit(Surv(time, status) ~ differentiation, data = df)
f8 <- survfit(Surv(time, status) ~ local_spread, data = df)
f9 <- survfit(Surv(time, status) ~ surg, data = df)
fits <- list(adhere = f1, rx = f2, sex =f3, obstruct = f4 , perfor=f5 , node4 =f6 , differentiation =f7 , local_spread =f8 , surg= f9)
# Visualize
legend.title <- list("adhere", "rx", "sex", "obstruct", "perfor", "node4", "differentiation", "local_spread", "surg")
ggsurvplot_list(fits, df, pval=TRUE) #legend.title = legend.title,$adhere
$rx
$sex
$obstruct
$perfor
$node4
$differentiation
$local_spread
$surg
attr(,"class")
[1] "list" "ggsurvplot_list"Figure 3: Survival curves stratified by categorical variables.
There is a statistically significant difference in survival times between two or more groups being compared as shown in the survival curves stratified by sex, obstruct, adherence to nearby organs, presence of four or more tumor positive nodes, differentiation, local spread and surgery. Survival time is not significantly different between male and female or among those with or without tumor perforation. The vertical lines on the Kaplan-Meier’s curve which show censoring data points are enriched after about 1,800 days, suggesting that many of the participants were lost to follow-up after this period. Patient characteristics that showed significant differences in survival time between two categories are important predictors of survival time. However, The Kaplan-Meier’s curve does not adjust for other variables, therefore, some of the predictors may not be significant when included in the multivariate Cox regression model.
IIId. Cox regression models
Model_1: Base Model
Model_1 <- coxph(Surv(time, status) ~ 1, data = df)
c_index_model_1 <- concordance(Model_1)
cat("Concordance of the base model:",c_index_model_1$concordance)Concordance of the base model: 0.5Model_2: Univariate analysis
# Univariate analysis
Model_2 <- coxph(Surv(time, status) ~ rx, data = df)
# Create the regression table and add concordance statistic
summary_table <- tbl_regression(Model_2, exponentiate = TRUE) %>%
add_glance_source_note(
label = list(concordance = "Concordance"),
include = c("concordance")
) %>%
modify_table_styling(
columns = p.value,
rows = p.value < 0.05,
text_format = "bold"
)
# Convert to gt table, increase font size, and adjust width
gt_table <- as_gt(summary_table) %>%
gt::tab_options(
table.font.size = px(17), # Increase font size
table.align = "center", # Align table to the left
table.width = pct(50) # Make the table wider
)
# Print the gt table
gt_table| Characteristic | HR1 | 95% CI1 | p-value |
|---|---|---|---|
| rx | |||
| Obs | — | — | |
| Lev | 0.97 | 0.78, 1.21 | 0.8 |
| Lev+5FU | 0.69 | 0.55, 0.87 | 0.002 |
| Concordance = 0.536 | |||
| 1 HR = Hazard Ratio, CI = Confidence Interval | |||
The coefficient of Lev is not significant, suggesting that there is no evidence that this treatment affects survival time compared to observation. however Lev+5Fu is significant (p=0.00175), indicating that the treatment Lev +5Fu has a statistically significant effect on survival time compared to the reference group. The negative sign indicates that this treatment group has a lower hazard and likely a longer survival time. The hazard ratio for Lex+5FU (0.690), indicating the risk of death is about 31% lower compared to the observation group.
Model_2: Cox proportional hazard assumption test
zph_test <- cox.zph(Model_2)
# Convert the Schoenfeld residuals test results to a data frame
zph_df <- as.data.frame(zph_test$table)
zph_df$Variable <- rownames(zph_df)
zph_df <- zph_df %>%
mutate(
chisq = round(chisq, 3),
p = round(p, 3)
)
zph_df %>%
kbl(caption = "Schoenfeld Residuals Test Results") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, position = "center") %>%
row_spec(0, bold = TRUE, align = "center")| chisq | df | p | Variable | |
|---|---|---|---|---|
| rx | 1.481 | 2 | 0.477 | rx |
| GLOBAL | 1.481 | 2 | 0.477 | GLOBAL |
# plot the Schoenfeld residuals
plot(zph_test)
The Schoenfeld residal plot shows that the residuals are scattered randomly and the smooth trend line is horizontal near 0. This suggests that the hazard ratio for rx (treatment status) is constant over time and the proportional hazard assumption is met. The global p-value is >0.05, indicating that the the assumption is met.
Model_3: All predictors
# Full variables: All predictors
Model_3<- coxph(Surv(time, status) ~ rx+ age + sex + perfor + adhere + surg + obstruct + differentiation + node4+ local_spread, data = df)
summary_table <- tbl_regression(Model_3, exponentiate = TRUE) %>%
add_glance_source_note(
label = list(concordance = "Concordance"),
include = c("concordance")
) %>%
modify_table_styling(
columns = p.value,
rows = p.value < 0.05,
text_format = "bold"
)
gt_table <- as_gt(summary_table) %>%
tab_options(
table.font.size = px(17), # Reduce font size
table.align = "center", # Align table to the center
table.width = pct(50), # Make the table wider
data_row.padding = px(2) # Reduce row padding
)
gt_table| Characteristic | HR1 | 95% CI1 | p-value |
|---|---|---|---|
| rx | |||
| Obs | — | — | |
| Lev | 0.98 | 0.79, 1.22 | 0.9 |
| Lev+5FU | 0.69 | 0.54, 0.87 | 0.002 |
| age | 1.01 | 1.00, 1.02 | 0.083 |
| sex | 1.04 | 0.86, 1.26 | 0.7 |
| perfor | 1.00 | 0.59, 1.70 | >0.9 |
| adhere | 1.18 | 0.92, 1.53 | 0.2 |
| surg | 1.27 | 1.03, 1.55 | 0.022 |
| obstruct | 1.33 | 1.06, 1.68 | 0.015 |
| differentiation | |||
| moderate | — | — | |
| poor | 1.43 | 1.13, 1.82 | 0.003 |
| well | 1.08 | 0.78, 1.50 | 0.6 |
| node4 | 2.55 | 2.10, 3.09 | <0.001 |
| local_spread | |||
| contiguous | — | — | |
| muscle | 0.39 | 0.23, 0.64 | <0.001 |
| serosa | 0.64 | 0.43, 0.94 | 0.023 |
| submucosa | 0.29 | 0.10, 0.83 | 0.021 |
| Concordance = 0.674 | |||
| 1 HR = Hazard Ratio, CI = Confidence Interval | |||
Model_3 summary shows that, when all variables are included in the model; age, sex, perforation and adherance are not significant predictors. wheras, rx, surg, obstruct, differentiation, node4 and local spread are significant predictors. The concordance of the multivariable model, 0.674, is higher than the univariate model (concordance =0.53), suggesting that the multivariate model is a better fit model. However, the model concordance did not improve when removing predictors that were not significant in Model_3.
Stepwise variable selection:
Stepwise variable selection approach was used to confirm that the significant predictors in Model_3 are the best set of predictors. The MASS package stepAIC() function was used for stepwise variable selection by using AIC (Akaike Information Criterion) as the measure to add or remove predictors from the model.
# Model_2 includes all variables
Model_2 <- coxph(Surv(time, status) ~ rx+ age + sex + perfor + adhere + surg + obstruct + differentiation + node4+
local_spread, data = df)
# stepwise selection
stepwise_model <- stepAIC(Model_2, direction = "both")Start: AIC=5741.4
Surv(time, status) ~ rx + age + sex + perfor + adhere + surg +
obstruct + differentiation + node4 + local_spread
Df AIC
- perfor 1 5739.4
- sex 1 5739.6
- adhere 1 5741.0
<none> 5741.4
- age 1 5742.5
- surg 1 5744.5
- obstruct 1 5745.1
- differentiation 2 5745.4
- rx 2 5749.9
- local_spread 3 5753.1
- node4 1 5822.0
Step: AIC=5739.4
Surv(time, status) ~ rx + age + sex + adhere + surg + obstruct +
differentiation + node4 + local_spread
Df AIC
- sex 1 5737.6
- adhere 1 5739.1
<none> 5739.4
- age 1 5740.5
+ perfor 1 5741.4
- surg 1 5742.5
- obstruct 1 5743.2
- differentiation 2 5743.5
- rx 2 5747.9
- local_spread 3 5751.1
- node4 1 5820.1
Step: AIC=5737.59
Surv(time, status) ~ rx + age + adhere + surg + obstruct + differentiation +
node4 + local_spread
Df AIC
- adhere 1 5737.3
<none> 5737.6
- age 1 5738.6
+ sex 1 5739.4
+ perfor 1 5739.6
- surg 1 5740.7
- obstruct 1 5741.2
- differentiation 2 5741.7
- rx 2 5746.2
- local_spread 3 5749.3
- node4 1 5818.2
Step: AIC=5737.26
Surv(time, status) ~ rx + age + surg + obstruct + differentiation +
node4 + local_spread
Df AIC
<none> 5737.3
+ adhere 1 5737.6
- age 1 5738.6
+ sex 1 5739.1
+ perfor 1 5739.2
- surg 1 5740.7
- obstruct 1 5740.9
- differentiation 2 5742.1
- rx 2 5746.0
- local_spread 3 5750.4
- node4 1 5817.4# Extract the coefficients of the selected model
selected_variables <- coef(stepwise_model)
# Print the selected variables
print(selected_variables) rxLev rxLev+5FU age
-0.010789186 -0.375746378 0.007365529
surg obstruct differentiationpoor
0.243870576 0.283164609 0.373782987
differentiationwell node4 local_spreadmuscle
0.069021620 0.929854405 -0.995556083
local_spreadserosa local_spreadsubmucosa
-0.501168855 -1.322018421 Based on a step-wise selection, rx, age, surg, obstruct, differentiation, node4 and local spread were selected for inclusion in the final model. Age is a significant predictor in the step-wise selected variables, but not in the full variable model (Model_2).
Model_4: Step-wise Selected variables
Model_4 <- coxph(Surv(time, status) ~ rx + age + surg + obstruct +
differentiation + node4 + local_spread, data = df)
summary_table <- tbl_regression(Model_4, exponentiate = TRUE) %>%
add_glance_source_note(
label = list(concordance = "Concordance"),
include = c("concordance")
) %>%
modify_table_styling(
columns = p.value,
rows = p.value < 0.05,
text_format = "bold"
)
gt_table <- as_gt(summary_table) %>%
tab_options(
table.font.size = px(17), # Reduce font size
table.align = "center", # Align table to the center
table.width = pct(50), # Make the table wider
data_row.padding = px(2) # Reduce row padding
)
gt_table| Characteristic | HR1 | 95% CI1 | p-value |
|---|---|---|---|
| rx | |||
| Obs | — | — | |
| Lev | 0.99 | 0.80, 1.23 | >0.9 |
| Lev+5FU | 0.69 | 0.54, 0.87 | 0.002 |
| age | 1.01 | 1.00, 1.02 | 0.069 |
| surg | 1.28 | 1.04, 1.56 | 0.018 |
| obstruct | 1.33 | 1.06, 1.67 | 0.015 |
| differentiation | |||
| moderate | — | — | |
| poor | 1.45 | 1.15, 1.84 | 0.002 |
| well | 1.07 | 0.77, 1.48 | 0.7 |
| node4 | 2.53 | 2.09, 3.07 | <0.001 |
| local_spread | |||
| contiguous | — | — | |
| muscle | 0.37 | 0.23, 0.61 | <0.001 |
| serosa | 0.61 | 0.41, 0.89 | 0.010 |
| submucosa | 0.27 | 0.09, 0.76 | 0.014 |
| Concordance = 0.672 | |||
| 1 HR = Hazard Ratio, CI = Confidence Interval | |||
The concordance of Model_4 did not improve when compared to Model_3, which included all predictors, but the model is less complex since in includes a smaller number of variables (see model comparison table below)
Model_4: Cox proportional hazard assumption test
# final model with step wise variable selection
zph_test <- cox.zph(Model_4)
# Convert the Schoenfeld residuals test results to a data frame
zph_df <- as.data.frame(zph_test$table)
zph_df$Variable <- rownames(zph_df)
zph_df <- zph_df %>%
mutate(
chisq = round(chisq, 3),
p = round(p, 3)
)
zph_df %>%
kbl(caption = "Schoenfeld Residuals Test Results") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, position = "center") %>%
row_spec(0, bold = TRUE, align = "center")| chisq | df | p | Variable | |
|---|---|---|---|---|
| rx | 2.335 | 2 | 0.311 | rx |
| age | 0.549 | 1 | 0.459 | age |
| surg | 0.020 | 1 | 0.888 | surg |
| obstruct | 6.148 | 1 | 0.013 | obstruct |
| differentiation | 17.459 | 2 | 0.000 | differentiation |
| node4 | 5.662 | 1 | 0.017 | node4 |
| local_spread | 7.083 | 3 | 0.069 | local_spread |
| GLOBAL | 37.525 | 11 | 0.000 | GLOBAL |
Model_4 does not meet Cox poportional hazard assumption as shown by p-values less than 0.05 in the Schoenfeld residuals test. Differentiation, node4 and obstruct variables did not meet proportional hazards assumption, suggesting that effects of these variables are not constant over time. To overcome this, stratification approach will be used, where the model will be stratified by variables that did not meet proportional hazard.
Model_5: Stratified Model
| Characteristic | HR1 | 95% CI1 | p-value |
|---|---|---|---|
| rx | |||
| Obs | — | — | |
| Lev | 0.98 | 0.79, 1.22 | 0.9 |
| Lev+5FU | 0.71 | 0.56, 0.89 | 0.003 |
| age | 1.01 | 1.00, 1.02 | 0.034 |
| surg | 1.30 | 1.06, 1.59 | 0.012 |
| node4 | 2.50 | 2.06, 3.04 | <0.001 |
| local_spread | |||
| contiguous | — | — | |
| muscle | 0.34 | 0.21, 0.56 | <0.001 |
| serosa | 0.58 | 0.39, 0.84 | 0.004 |
| submucosa | 0.24 | 0.08, 0.67 | 0.007 |
| Concordance = 0.674 | |||
| 1 HR = Hazard Ratio, CI = Confidence Interval | |||
Model_5 Cox proportional hazard assumption test
zph_test <- cox.zph(Model_5)
zph_test chisq df p
rx 2.0007 2 0.368
age 0.6704 1 0.413
surg 0.0142 1 0.905
node4 4.2882 1 0.038
local_spread 5.2976 3 0.151
GLOBAL 12.4113 8 0.134# Convert the Schoenfeld residuals test results to a data frame
zph_df <- as.data.frame(zph_test$table)
zph_df$Variable <- rownames(zph_df)
zph_df <- zph_df %>%
mutate(
chisq = round(chisq, 3),
p = round(p, 3)
)
zph_df %>%
kbl(caption = "Schoenfeld Residuals Test Results") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, position = "center") %>%
row_spec(0, bold = TRUE, align = "center")| chisq | df | p | Variable | |
|---|---|---|---|---|
| rx | 2.001 | 2 | 0.368 | rx |
| age | 0.670 | 1 | 0.413 | age |
| surg | 0.014 | 1 | 0.905 | surg |
| node4 | 4.288 | 1 | 0.038 | node4 |
| local_spread | 5.298 | 3 | 0.151 | local_spread |
| GLOBAL | 12.411 | 8 | 0.134 | GLOBAL |
# plot the Schoenfeld residuals
plot(zph_test)
After model stratification by obstruct and differentiation, the proportional hazard assumption is met as the global p-value is >0.05. Node4 slightly violates assumption, but the final model is not stratified by node4 because the model concordance is attenuated when stratifying by node4.
Evaluate multicollinearity of variables in Model_5
vif <- vif(Model_5)
print(vif) GVIF Df GVIF^(1/(2*Df))
rx 1.008066 2 1.002010
age 1.019142 1 1.009525
surg 1.008430 1 1.004206
strata(obstruct) 3.071784 0 Inf
strata(differentiation) 3.071784 0 Inf
node4 1.023390 1 1.011628
local_spread 1.017073 3 1.002825None of the variables has VIF >5, therefore multicollinearity is not a problem in Model_5
Model comparison
# Fit the models and store them in a list
models <- list(Model_1, Model_2, Model_3, Model_4, Model_5)
# Add descriptions for each model
descriptions <- c(
"Base model",
"Treatment",
"Full variables",
"Stepwise-selected variables",
"Stratified"
)
# Create a data frame to store results
results <- data.frame(
Model = character(),
Description = character(),
AIC = numeric(),
BIC = numeric(),
C_Index = numeric(),
stringsAsFactors = FALSE
)
# Function to calculate and store metrics for each model
for (i in seq_along(models)) {
model <- models[[i]]
# Extract AIC and BIC
aic <- AIC(model)
bic <- BIC(model)
# Add C-index
c_index <- concordance(model)$concordance
# Append results to the data frame
results <- rbind(results, data.frame(
Model = paste("Model", i),
Description = descriptions[i],
AIC = aic,
BIC = bic,
C_Index = round(c_index, 3)
))
}
# Print the table using kable and kableExtra
results %>%
kbl(caption = "Model Evaluation Metrics") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, position = "center") %>%
column_spec(1, bold = TRUE) %>%
column_spec(2:5, width = "10em") %>%
kable_styling(position = "center")| Model | Description | AIC | BIC | C_Index |
|---|---|---|---|---|
| Model 1 | Base model | 5860.383 | 5860.383 | 0.500 |
| Model 2 | Treatment | 5741.401 | 5798.993 | 0.674 |
| Model 3 | Full variables | 5741.401 | 5798.993 | 0.674 |
| Model 4 | Stepwise-selected variables | 5737.261 | 5782.511 | 0.672 |
| Model 5 | Stratified | 4567.829 | 4600.739 | 0.674 |
Based on the model evaluation metrics, Model_5 has the smallest AIC and BIC values, and achieved highest concordance (c-index). Although Model_2 and Model_5 have the same concordance, Model_5 is the best model as indicated by the smaller AIC and BIC values when compared to Model_2. 5-fold cross-validation was performed to test the robustness of Model_5.
K-fold cross validation
set.seed(1234)
# Cox model
cox_model <- coxph(Surv(time, status) ~ rx + age + surg + strata(obstruct) + strata(differentiation) + node4 +
local_spread, data = df)
# calculate the original c-index
c_index_original <- concordance(cox_model)
cat("original c-index:", c_index_original$concordance, "\n")original c-index: 0.674316 # create a function for calculating c-index in each fold using concordance()
cox_cindex <- function(train_data, test_data) {
fit <- coxph(Surv(time, status) ~ rx + age + surg + strata(obstruct) + strata(differentiation) + node4 + local_spread, data = train_data)
# Calculate concordance on test data
c_index <- concordance(fit, newdata = test_data)$concordance
return(c_index)
}
# perform 5-fold cross-validation with stratification
K <- 5
folds <- createFolds(c(df$status, df$differentiation, df$rx), k = K, list = TRUE, returnTrain = TRUE)
cv_c_indices <- sapply(folds, function(train_indices) {
train_data <- df[train_indices, ]
test_data <- df[-train_indices, ]
cox_cindex(train_data, test_data) # use the concordance function inside cox_cindex
})
# Print cross-validated c-indices
print(cv_c_indices) Fold1 Fold2 Fold3 Fold4 Fold5
0.7156051 0.6605045 0.6562016 0.6477106 0.6943820 # cross-validation c-indices
cat("cross-validated c-Indices for each fold:", cv_c_indices, "\n")cross-validated c-Indices for each fold: 0.7156051 0.6605045 0.6562016 0.6477106 0.694382 cat("mean cross-validated c-Index:", mean(cv_c_indices), "\n")mean cross-validated c-Index: 0.6748808 # plot cross-validation c-indices
plot(cv_c_indices, type = "b", xlab = "Fold", ylab = "c-index", main = "c-index across folds")
Model_5 c-index (0.674) and mean cross-validation c-index (0.675) is very similar, suggesting the the final stratified model is stable and is not overfitting.
IV. Conclusions
- Treatment with Levamisole + 5-FU decreases the hazard of death from colon cancer by 29.5% (HR =0.71, 95% CI: 0.56-0.89; p=0.0035).
- Having more than 4 tumor positive lymph nodes significantly increases the hazard of death by 150.3% (p<0.0001).
- Having a long wait period from surgery to registration for trial is associated with an increase in the hazard by 99.5% (p=0.01).
- Patients with local tumor spread to the submucosa, muscle and serosa have a reduction in the hazard by 76% (p=0.0007), 66% (p<0.001), and 43% (0.004), respectively, compared to those with contiguous organ spread.
- The concordance of the model (0.67) indicates moderate predictive accuracy for survival time.
- Other variables that were not included in the study may contribute to survival time.
References
[1]
S. J. Walters, “Analyzing time to event outcomes with a cox regression model,” Wiley Interdiscip. Rev. Comput. Stat., vol. 4, no. 3, pp. 310–315, May 2012.
[2]
V. Patil and S. Dessai, “Testing and interpreting assumptions of COX regression analysis,” Cancer Res. Stat. Treat., vol. 2, no. 1, p. 108, 2019.
[3]
V. Backmann et al., “Comprehensive strain assessment and mortality after acute myocardial infarction: A retrospective observational study based on the essen coronary artery disease registry,” Heart, Sep. 2024.
[4]
Z. Zhang, J. Reinikainen, K. A. Adeleke, M. E. Pieterse, and C. G. M. Groothuis-Oudshoorn, “Time-varying covariates and coefficients in cox regression models,” Ann. Transl. Med., vol. 6, no. 7, pp. 121–121, Apr. 2018.
[5]
T. Therneau and P. Grambsch, Modeling survival data: Extending the cox model, 1st ed. in Statistics for biology and health. New York, NY: Springer, 2010.
[6]
D. G. Kleinbaum, “The cox proportional hazards model and its characteristics,” in Survival analysis, New York, NY: Springer New York, 1996, pp. 83–128.
[7]
T. Intrieri, G. Manneschi, and A. Caldarella, “10-year survival in female breast cancer patients according to ER, PR and HER2 expression: A cancer registry population-based analysis,” J. Cancer Res. Clin. Oncol., vol. 149, no. 8, pp. 4489–4496, Jul. 2023.
[8]
Y. Wu et al., “Effect of adjuvant chemotherapy on the survival outcomes of elderly breast cancer: A retrospective cohort study based on SEER database,” J. Evid. Based Med., vol. 15, no. 4, pp. 354–364, Dec. 2022.
[9]
A. Ptak-Chmielewska and J. P. E. Gonzalez, “Default prediction using the cox regression model and macroeconomic conditions – a lifetime perspective,” Econometrics, vol. 28, no. 2, pp. 50–61, 2024.
[10]
K. K.-K. Wong, “Using cox regression to model customer time to churn in the wireless telecommunications industry,” J. Target. Meas. Anal. Mark., vol. 19, no. 1, pp. 37–43, Mar. 2011.
[11]
M. I. Rodrı́guez-Borbón, M. A. Rodrı́guez-Medina, L. A. Rodrı́guez-Picón, A. Alvarado-Iniesta, and N. Sha, “Reliability estimation for accelerated life tests based on a cox proportional hazard model with error effect,” Qual. Reliab. Eng. Int., vol. 33, no. 7, pp. 1407–1416, Nov. 2017.
[12]
L. Cheng, W. K. Chan, L. Zhu, M. H. Chao, and Y. Wang, “Confronting inequalities and bridging the divide: A retrospective study assessment of country-level COVID-19 vaccine equality with a cox regression model,” Vaccines (Basel), vol. 12, no. 5, p. 552, May 2024.
[13]
T. M. Therneau, A package for survival analysis in r. 2024. Available: https://CRAN.R-project.org/package=survival
[14]
Terry M. Therneau and Patricia M. Grambsch, Modeling survival data: Extending the Cox model. New York: Springer, 2000.
[15]
C. G. Moertel et al., “Levamisole and fluorouracil for adjuvant therapy of resected colon carcinoma,” N. Engl. J. Med., vol. 322, no. 6, pp. 352–358, Feb. 1990.
[16]
A. Kassambara, M. Kosinski, and P. Biecek, Survminer: Drawing survival curves using ’ggplot2’. 2021. Available: https://CRAN.R-project.org/package=survminer
[17]
W. N. Venables and B. D. Ripley, Modern applied statistics with s, Fourth. New York: Springer, 2002. Available: https://www.stats.ox.ac.uk/pub/MASS4/
[18]
R Core Team, R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing, 2023. Available: https://www.R-project.org/
[19]
A. C. Davison and D. V. Hinkley, Bootstrap methods and their applications. Cambridge: Cambridge University Press, 1997. Available: doi:10.1017/CBO9780511802843
[20]
Angelo Canty and B. D. Ripley, Boot: Bootstrap r (s-plus) functions. 2024.